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In this work we reexamine the LDA+U method of Anisimov and coworkers in the framework of 
a plane-wave pseudopotential approach. A simplified rotational-invariant formulation is adopted. 
The calculation of the Hubbard U entering the expression of the functional is discussed and a 
linear response approach is proposed that is internally consistent with the chosen definition for 
the occupation matrix of the relevant localized orbitals. In this way we obtain a scheme whose 
functionality should not depend strongly on the particular implementation of the model in ab-initio 
calculations. We demonstrate the accuracy of the method, computing structural and electronic 
properties of a few systems including transition and rare-earth correlated metals, transition metal 
monoxides and iron-silicate. 



INTRODUCTION 

The description and understanding of electronic prop- 
erties of strongly correlated materials is a very impor- 
tant and long standing problem for ab-initio calculations. 
Widely used approximations for the exchange and corre- 
lation energy in density functional theory (DFT), mainly 
based on parametrization of (nearly) homogeneous elec- 
tron gas, miss important features of their physical behav- 
ior. For instance both local spin-density approximation 
(LSDA) and spin-polarized generalized gradient approx- 
imations (a— GGA), in their several flavors, fail in pre- 
dicting the insulating behavior of many simple transition 
metal oxides (TMO), not only by severely underestimat- 
ing their electronic band gap but, in most cases, produc- 
ing a qualitatively wrong metallic ground state. 

TMOs have represented for long time the most notable 
failure of DFT. When the high-T c superconductors en- 
tered the scene (their parent materials are also strongly 
correlated systems) the quest for new approaches that 
could describe accurately these systems by first princi- 
ples received new impulse, and in the last fifteen years 
many methods were proposed in this direction. Among 
these, LDA+U approach, first introduced by Anisimov 
and coworkers has allowed to study a large vari- 

ety of strongly correlated compounds with considerable 
improvement with respect to LSDA or a— GGA results. 
The successes of the method have led to further develop- 
ments during the last decade which have produced very 
sophisticated theoretical approaches |^ and efficient nu- 
merical techniques. 

The formal expression of LDA+U energy functional 
is adapted from model hamiltonians (Hubbard model 
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in particular) that represent the "natural" theoretical 
framework to deal with strongly correlated materials. As 
in these models, a small number of localized orbitals is 
selected and the electronic correlation associated to them 
is treated in a special way. The obtained results strongly 
depend on the definition of the localized orbitals and on 
the choice of the interaction parameters used in the cal- 
culation, that should be determined in an internally con- 
sistent with. This is not always done and a widespread 
but, in our opinion, unsatisfactory attitude is to deter- 
mine the value of the electronic couplings by seeking a 
good agreement of the calculated properties with the ex- 
perimental results in a semiempirical way. 

In this work a critical reexamination of the LDA+U 
approach is proposed, which starting from the formula- 
tion of Anisimov and coworkers 0, @] > an d its further 
improvements 0, 0> 0> develops a simpler approxima- 
tion. This is, in our opinion, the "minimal" extension 
of the usual approximate DFT (LDA or GGA) schemes 
needed when atomic-like features are persistent in the 
solid environment. 

In the central part of this work we describe a method, 
based on a linear response approach, to calculate in 
an internally consistent way — without aprioristic as- 
sumption about screening and/or basis set employed in 
the calculation — the interaction parameters entering the 
LDA+U functional used. In this context our plane-wave 
pseudopotential (PWPP) implementation of the LDA+U 
approach is presented and discussed in some details. We 
stress however that the proposed method is basis-set in- 
dependent. 

Our methodology is then applied to the study of the 
electronic properties of some real materials, chosen as 
representative of "normal" (bulk iron) and correlated 
(bulk cerium) metals, as well as a few examples of 
strongly correlated systems (iron oxide, nickel oxide and 
fayalite) . 
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STANDARD LDA+U IMPLEMENTATION: 

In order to account explicitly for the on-site Coulomb 
interaction responsible for the correlation gap in Mott in- 
sulators and not treated faithfully within LDA, Anisimov 
and coworkers |il LJ l| correct the standard functional 
adding an on-site Hubbard-like interaction, Enub- 
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where n(r) is the electronic density, and are the 
atomic-orbital occupations for the atom / experiencing 
the "Hubbard" term. The last term in the above equa- 
tion is then subtracted in order to avoid double count- 
ing of the interactions contained both in En u b and, in 
some average way, in Eld a- In this term the total, spin- 
projected, occupation of the localized manifold is used: 

In its original definition the functional defined in Eq. 
^was not invariant under rotation of the atomic-orbital 
basis set used to define the occupancies n^- A rotation- 
ally invariant formulation has then been introduced 0, |(J 
where the orbital dependence of En u b is borrowed from 
atomic Hartree-Fock with renormalized slater integrals: 



by assuming atomic values for F A /F 2 and F G /F 4 ratios. 

To obtain U and J, Anisimov and coworkers 0,0 pro- 
pose to perform LMTO calculations in supercells in which 
the occupation of the localized orbitals of one atom is 
constrained. The localized orbitals of all atoms in the 
supercell are decoupled from the remainder of the basis 
set. This makes the treatment of the local orbitals an 
atomic-like problem — making it easy to fix their occupa- 
tion numbers — and allows to use Janak theorem [9| to 
identify the shift in the corresponding eigenvalue with 
the second-order derivative of the LDA total energy with 
respect to orbital occupation. It has however the effect of 
leaving a rather artificial system to perform the screen- 
ing, in particular when it is not completely intra-atomic. 
In elemental metallic Iron, for instance, Anisimov and 
Gunnarsson y| showed that only half of the screening 
charge is contained in the Wigner-Seitz cell. This fact, 
in addition to a sizable error due to the Atomic Sphere 
Approximation used @ , could be at the origin of the se- 
vere overestimation of the computed on-site coulomb in- 
teraction with respect to estimates based on comparison 
of spectroscopic data and model calculations [1 (A 1 1 lj . 



EHub[W mm ,}} 



E 

{m},cr,7 



{(' 



-{(m,m"\V e , 
-(to, m"\V ee 



m 



\V ee \n 

,to"') 
, m'))n 



')n Ia 

/"■mrr 



BASIS SET INDEPENDENT FORMULATION OF 
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where I is the angular moment of the localized (d or /) 
electrons and 
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The double-counting term E c i c is given by: 
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The radial Slater integrals F k are the parameters of the 
model (F°,F 2 and F A for d electrons, while also F 6 must 
be specified for / states) and are usually re-expressed 
in terms of only two parameters, U and J, describing 
screened on-site Coulomb and exchange interaction, 
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Some aspects of currently used LDA+U formulation, 
and in particular of the determination of the parameters 
entering the model, have been so far tied to the LMTO 
approach. This is not a very pleasant situation and some 
efforts have been done recently 0, E3 to reformulate the 
method for different basis sets. Here we want to elabo- 
rate further on these attempts and provide an internally 
consistent, basis-set independent, method for the calcu- 
lation of the needed parameters. 



Localized orbital occupations 

In order to fully define how the approach works the first 
thing to do is to select the degrees of freedom on which 
"Hubbard W will operate and define the corresponding 
occupation matrix, n^ m j . Although it is usually straight- 
forward to identify in a given system the atomic levels to 
be treated in a special way (the d electrons in transition 
metals and the / ones in the rare earths and actinides 
series) there is no unique or rigorous way to define occu- 
pation of localized atomic levels in a multi-atom system. 
Equally legitimate choices for n^ m , are i) projections on 
normalized atomic orbitals, or ii) projections on Wannier 
functions whenever the relevant orbitals give raise to iso- 
lated band manifolds, or Hi) Mulliken population or iv) 
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integrated values in (spherical) regions around the atoms 
of the angular-momentum-decomposed charge densities. 
Taking into account the arbitrariness in the definition of 
n mm' n0 particular significance should be attached to any 
of them (or other that could be introduced) and the use- 
fulness and reliability of an approximate DFT+U method 
(aDFT+U), and of its more recent and involved evolu- 
tions like the aDFT+DMFT method, should be judged 
from its ability to provide a correct physical picture of 
the systems under study irrespective of the details of the 
formulation, once all ingredients entering the calculation 
are determined consistently. 

All above mentioned definitions for the occupation ma- 
trices can be put in the generic form 



i' — /ku(V'kJ-Pmm'IV'kij) 



(5) 
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where tp£ v ^ s ^ ne va lence electronic wavefunction corre- 
sponding to the state (ku) with spin a of the system and 
/£„ is the corresponding occupation number. The P^m' 's 
are generalized projection operators on the localized- 
electron manifold that satisfy the following properties: 
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In particular P 1 = J2 m Pmm ^ s the projector on the com- 
plete manifold of localized states associated with atom at 
site I and therefore 
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is the total localized-states occupation for site I. Orthog- 
onality of projectors on different sites is not assumed. 

In the applications discussed in this work we will define 
localized-level occupation matrices projecting on atomic 
pseudo-wavefunctions. The needed projector operators 
are therefore simply 



(8) 



where is the valence atomic orbital with angular 

momentum component \lm) of the atom sitting at site I 
(the same wavefunctions are used for both spins). Since 
we will be using ultrasoft pseudopotentials to describe 
valence-core interaction, all scalar products between crys- 
tal and atomic pseudo-wavefunctions are intended to in- 
clude the usual S matrix describing orthogonality in pres- 
ence of charge augmentation [lflj . 

As already mentioned, other choices could be used as 
well and different definitions for the occupation matrices 
will require, in general, different values of the parameter 
entering the aDFT+U functional, as it has been pointed 
out recently also by Pickett et al. Q where, for instance, 
the value of Hubbard U in FeO shifts from 4.6 to 7.8 eV 
when atomic d-orbitals for Fe 2+ ionic configuration are 



used instead of those of the neutral atom. In an early 
study 19] the U parameter in La2Cu04 varies from 6.8 
to 7.7 eV upon variation of the atomic sphere radius em- 
ployed in the LMTO calculation. As pointed out in these 
works it is not fruitful to compare numerical values of 
U obtained by different methods but rather comparison 
should be made between results of complete calculations. 



A simplified rotationally invariant scheme and the 
meaning of U 

In order to simplify our analysis and gaining a more 
transparent physical interpretation of the " +U" correc- 
tion to standard aDFT functionals we concentrate on the 
main effect associated to on-site Coulomb repulsion. We 
thus neglect the important but somehow secondary ef- 
fects associated to non sphericity of the electronic inter- 
action and the proper treatment of magnetic interaction, 
that in the currently used rotational invariant method is 
dealt with assuming a screened Hartree-Fock form. |5| . 

We are therefore going to assume in the following that 
parameter J describing these effects can be set to zero, or 
alternatively that its effects can be mimicked redefining 
the U parameter as U e ff = U — J, a practice that have 
been sometime used in the literature 

The Hubbard correction to the energy functional, Eqs. 
13 and greatly simplifies and reads: 
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Choosing for the localized orbitals the representation 
that diagonalizes the occupation matrices 



(10) 



with < < 1, the energy correction becomes 

Eu[{n!Z m ,}\ = lEE^ 1 -^)- ( H ) 

I,i7 i 

from where it appears clearly that the energy correction 
introduces a penalty, tuned by the value of the U param- 
eter, for partial occupation of the localized orbitals and 
thus favors disproportionation in fully occupied (A ~ 1) 
or completely empty (A rs 0) orbitals. This is the ba- 
sic physical effect built in the aDFT+U functional and 
its meaning can be traced back to known deficiencies of 
LDA or GGA for atomic systems. 

An atom in contact with a reservoir of electrons can 
exchange integer numbers of particles with its environ- 
ment. The intermediate situation with fractional number 



4 



E(N+2) 



<D 

c 

<D 

-4— » 

O 




N N+1 

Number of electrons 

FIG. 1: Sketch of the total energy profile as a function of 
number of electrons in a generic atomic system in contact 
with a reservoir. The bottom curve is simply the difference 
between the other two (the LDA energy and the "exact" result 
for an open system). 



of electrons in this open atomic system is described not 
by a pure state wave function, but rather by a statisti- 
cal mixture so that, for instance, the total energy of a 
system with N + u> electrons (where N is an integer and 
< u> < 1) is given by: 
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where En and En+i are the energies of the system cor- 
responding to states with N and N+1 particles respec- 
tively, while u) represents the statistical weight of the 
state with N+1 electrons. The total energy of this open 
atomic system is thus represented by a series of straight- 
line segments joining states corresponding to integer oc- 
cupations of the atomic orbitals as depicted in fig. ^ The 
slope of the energy vs electron-number curve is instead 
piece-wise constant, with discontinuity for integer num- 
ber of electrons, and corresponds to the electron affinity 
(ionization potential) of the N (N+1) electron system. 

Exact DFT correctly reproduce this behavior 0, ^| , 
which is instead not well described by the LDA or GGA 
approach, which produces total energy with unphysical 
curvatures for non integer occupation and spurious min- 
ima in correspondence of fractional occupation of the or- 
bital of the atomic system. This leads to serious problems 
when one consider the dissociation limit of hetero-polar 
molecules or an open-shell atom in front of a metallic sur- 
face an d is at the heart of the LDA/GGA failure 
in the description of strongly correlated systems The 
unphysical curvature is associated basically to the incor- 
rect treatment by LDA or GGA of the self-interaction 
of the partially occupied Kohn-Sham orbital that gives a 



non-linear contribution to the total energy with respect 
to orbital occupation (with mainly a quadratic term com- 
ing from the Hartree energy not canceled properly in the 
exchange-correlation term). 

Nevertheless, it is well known that total energy dif- 
ferences between different states can be reproduced quite 
accurately by the LDA (or GGA) approach, if the oc- 
cupation of the orbitals is constrained to assume integer 
values. As an alternative, we can recover the physical sit- 
uation (an approximately piece-wise linear total energy 
curve) by adding a correction to the LDA total energy 
which vanishes for integer number of electrons and elim- 
inates the curvature of the LDA energy profile in every 
interval with fractional occupation (bottom curve of fig. 
01 . But this is exactly the kind of correction that is pro- 
vided by eq. El if the numerical value of the parameter U 
is set equal to the curvature of the LDA (GGA) energy 
profile. 

This clarifies the meaning of the interaction parame- 
ter U as the (unphysical) curvature of the LDA energy 
as a function of N which is associated with the spuri- 
ous self-interaction of the fractional electron injected into 
the system. From this analysis it is clear that the nu- 
merical value of U will depend in general not only, as 
noted in the preceding section, on the definition adopted 
for the occupation matrices but also on the particular 
approximate exchange-correlation functional to be cor- 
rected, and should vanish if the exact DFT functional 
were used. 

The situation is of course more complicated in solids 
where fractional occupations of the atomic orbitals can 
occur due to hybridization of the localized atomic-like 
orbitals with the crystal environment and the unphysi- 
cal part of the curvature has to be extracted from the 
total LDA/GGA energy, which contains also hybridiza- 
tion effects. In the next section this problem is discussed 
and a linear response approach to evaluate Hubbard U 
is proposed. 



Internally consistent calculation of U 



Following previous seminal works [3, LLa, LUj we com- 
pute U by means of constrained-density-functional cal- 
culations [jof. What we need is the total energy as a 
function of the localized-level occupations of the "Hub- 
bard" sites: 



E[{ qi }}= mm \E[n(r)} 



^ar(nr - qi) 
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where the constraints on the site occupations, n/'s from 
Eq. are applied employing the Lagrange multipliers, 
a/'s. From this dependence we can compute numeri- 
cally the curvature of the total energy with respect to 
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the variation, around the unconstrained values {rij} , 
of the occupation of one isolated site. A supercell ap- 
proach is adopted in which occupation of one representa- 
tive site in a sufficiently large supercell is changed leav- 
ing unchanged all other site occupations. This curvature 
contains the energy cost associated to the localization of 
an electron on the chosen site including all screening ef- 
fects from the crystal environment, but it is not yet the 
Hubbard U we want to compute. In fact, had we com- 
puted the same quantity from the total energy of the non- 
interacting Kohn-Sham problem associated to the same 
system, 



E KS [{qi}}= min {E KS [n{r)}+ V«f S ( 

n(r),a, I 



rii - qi) 



(14) 

we would have obtained a non vanishing results as well 
because by varying the site occupation a rehybridization 
of the localized orbitals with the other degrees of freedom 
is induced that gives rise to a non-linear change in the 
energy of the system. This curvature coming from re- 
hybridization, originating from the non-interacting band 
structure but present also in the interacting case, has 
clearly nothing to do with the Hubbard U of the inter- 
acting system and should be subtracted from the total 
curvature: 



U = 



d 2 E[{ qi }] d 2 E KS [{ qi }) 



dq] 



dqj 
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In Ref. y| Anisimov and Gunnarsson, in order to avoid 
dealing with the above mentioned non-interacting curva- 
ture, exploited the peculiarities of the LMTO method, 
used in their calculation, and decoupled the chosen lo- 
calized orbitals from the remainder of the crystal by sup- 
pressing in the LMTO hamiltonian the corresponding 
hopping terms. This reduced the problem to the one of 
an isolated atom embedded in an artificially disconnected 
charge background. Thanks to Janak theorem [9j the 
second order derivative of the total energy in Eg. 1 131 can 
then be recast as a first order derivative of the localized- 
level eigenvalue. In our approach the role played in Refs. 
0> 13 by the eigenvalue of the artificially isolated atom 
is taken by the Lagrange multiplier, used to enforce level 
occupation 20] : 
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daf s 
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At variance with the original method of Refs. m 
our approach we need to compute and subtract the band- 
structure contribution, —daf s /dqi, from the total cur- 
vature, but, in return, Hubbard U is computed in exactly 
the same system to which it is going to be applied and 



the screening from the environment is more realistically 
included. The present method was inspired by the linear 
response scheme proposed by Pickett and coworkers Q 
where however the role of the non-interacting curvature 
was not appreciated. 

In actual calculations constraining the localized orbital 
occupations is not very practical and it is easier to pass, 
via a Legendre transform, to a representation where the 
independent variables are the ar's 
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Variation of these functionals with respect to wavefunc- 
tions shows that the effect of the a/'s is to add to the 
single particle potential a term, AV = J2i ctiP 1 (or 
Ay = J^i a i P 1 f° r the non-interacting case), where 
localized potential shifts of strength ai (af s ) are ap- 
plied to the localized levels associated to site /. 

It is useful to introduce the (interacting and non- 
interacting) density response functions of the system with 
respect to these localized perturbations: 
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Using this response-function language, the effective in- 
teraction parameter U associated to site / can be recast 
as: 



U 
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(19) 



that is reminiscent of the well known random-phase ap- 
proximation |2l) in linear response theory giving the in- 
teracting density response in terms of the non-interacting 
one and the Coulomb kernel. A similar result is obtained 
within DFT linear response [13] where the interaction 
kernel also contains an exchange-correlation part. 

The response functions, Eq. El needed in Eq. El 
are computed taking numerical derivatives. We per- 
form a well converged LDA calculation for the uncon- 
strained system (ai = for all sites in the supercell) 
and — starting from its self-consistent potential — we add 
small (positive and negative) potential shifts on each 
non equivalent "Hubbard" site J and compute the vari- 
ation of the occupations, n.j's, for all sites in the super- 
cell in two ways: i) letting the Kohn-Sham potential of 
the system readjust self-consistently to optimally screen 
the localized perturbation, Ay = ajPj, and ii) with- 
out allowing this screening. This latter result is nothing 
but the variation computed from the first iteration in 
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the self-consistent cycle leading eventually to the former 
(screened) results. The site-occupation derivatives calcu- 
lated according to i) and ii) give the matrices xu an d 
X?j respectively. 

Further considerations 

Before moving to examine some specific examples in 
the next section, let's end the present one discussing a 
few additional technical points. 

As mentioned earlier, Hubbard U is computed, ideally, 
from variation of the site occupation of a single site in an 
infinite crystal and in practice adopting a supercell ap- 
proach where periodically repeated sites are perturbed 
coherently. In order to speed up the convergence of the 
computed U with supercell size it may result useful to 
enforce explicitly charge neutrality for the perturbation, 
that is to be introduced in the response functions, thus 
enhancing its local character and reduce the interaction 
with its periodic images. In this procedure we introduce 
in the response functions, x an d X°j — m addition to 
the degrees of freedom associated to the localized sites — 
also a "delocalized background" representing all other 
degrees of freedom in the system. This translates in one 
more column and row in the response matrices, whose el- 
ements are determined imposing overall charge neutrality 
of the perturbed system for all localized perturbations, 
(E/X/J = °, E/Xjj = 0, VJ) and absence of any 
charge density variation upon perturbing the system with 
a constant potential (J2j XlJ = 0, J2j X?j = 0) V/). 
From a mathematical point of view both x an d Xo acquire 
a null eigenvalue, corresponding to a constant potential 
shift, and the needed inversions in Eg. 1191 must be taken 
with care. It can be shown that their singularities cancel 
out when computing the difference Xq 1 — X 1 an d the 
final result is well defined. We stress that in the limit 
of infinitely large supercell the coupling with the back- 
ground gives no contribution to the computed U, but we 
found that this limit is approached more rapidly when 
this additional degrees of freedom is included. 

In the same spirit we found that the spatial locality 
of the response matrices can be rather different from the 
one of their inverse and a supercell sufficient to decou- 
ple the periodically repeated response may be too small 
to describe correctly the inverse in eq. As a prac- 
tical procedure, therefore, after evaluating the response 
function matrices in a given supercell, we extrapolate the 
result to much larger supercells assuming that the most 
important matrix elements in xa an d X involve the atoms 
in the few nearest coordination-shells accessible in the 
original supercell. The corresponding matrix elements of 
the larger supercell are filled with the values extracted 
from the smaller one while all other, more distant, inter- 
actions are neglected. Again, when a sufficiently large 
supercell to extract the matrix elements of the response 



functions is considered, the effect of this extrapolation 
vanishes, but, as we will see in the following, this scheme 
capture a large fraction of the system-size dependence of 
the calculated U and it may allow to reach more rapidly 
the converged result. 

As a final remark we notice that the electronic struc- 
ture of a system described within the LDA+U approach 
may largely differ from the one obtained within the LDA 
used to compute U. In a more refined approach one 
might seek internal consistency between the band struc- 
ture used in the calculation of U and the one obtained 
using it. We have not addressed this issue here, but one 
can imagine performing the same type of analysis leading 
to the U determination for a functional already contain- 
ing an LDA+U correction. The computed U would in 
that case be a correction to be added to the original U 
and internal consistency would be reached when the cor- 
rection vanishes. 

EXAMPLES 
Metals: Iron and Cerium 

In their seminal paper Anisimov and Gunnarsson || 
computed the effective on site Coulomb interaction be- 
tween the localized electrons in metallic Fe and Ce. For 
Ce the calculated Coulomb interaction was about 6 eV 
in good agreement with empirical and experimental esti- 
mates ranging from 5 to 7 eV HiJ, 0] , while the result 
for Fe (also about 6 eV) was surprisingly high since U was 
expected to be in the range of 1-2 eV for elemental tran- 
sition metals, with the exception of Ni Let us 
apply the present approach to these two system, starting 
with Iron. 

In its ground state elemental Iron has a ferromag- 
netic (FM) spin arrangement and a body-centered cu- 
bic (BCC) structure. Gradient corrected exchange- 
correlation functional are needed in order to stabilize the 
experimental structure as compared with non-magnetic 
face-centered cubic (FCC) structure preferred by LDA. 
The Perdew-Burke-Ernzherof (PBE) ^ GGA func- 
tional was employed here. Iron ions were represented 
by ultrasoft pseudopotential and kinetic energy cutoffs 
of 35 Ry and 420 Ry were adopted for wavefunction and 
charge density Fourier expansion. Brillouin Zone inte- 
grations where performed using 8x8x8 Monkhorst and 
Pack special point grids (2(| using Methfessel and Paxton 
smearing technique [23] with a smearing width of 0.005 
Ry in order to smooth the Fermi distribution. 

The calculation of the effective Hubbard U followed the 
procedure outlined in preceding section: a supercell was 
selected containing a number of inequivalent Iron atoms; 
then, after a well converged self-consistent calculation, 
we applied to one of these atoms small, positive and neg- 
ative, potential shifts, AV = aPd (with a — ±0.2-0.5 
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FIG. 2: Calculated Hubbard U in metallic Iron for different 
supercells. Lines connect results from the cell-extrapolation 
procedure described in the text and different symbols corre- 
spond to inclusion of screening contributions up to the indi- 
cated shell of neighbors of the perturbed atom. 



eV), where is the projector on the localized d electron 
of the selected atom. From the variation of the d-level 
occupations of all Iron atoms in the cell one column of 
X and x° response functions was extracted and all other 
matrix elements were reconstructed by symmetry, includ- 
ing the background as explained previously. Hubbard U 
was then calculated from Eq. HT?1 

In order to describe response for an isolated pertur- 
bation four supercells were considered: i) a simple cu- 
bic (SC) cell containing two inequivalent iron atoms, the 
perturbed atom and one of its nearest neighbors; ii) 
a 2x2x2 BCC supercell containing 8 inequivalent Iron 
atoms, 4 in the nearest-neighbor shell of the perturbed 
atom and 3 belonging to the second shell of neighbors; 
Hi) a 2x2x2 SC cell containing 16 atoms, including also 
some third nearest-neighbor atom and iv) a4x4x4 BBC 
supercell containing 64 inequivalent Iron atoms; we used 
this largest cell just to extrapolate the results from the 
smaller ones. 

The convergence properties of the effective U of bulk 
iron with the size of the used supercell are shown in fig.|2 

The Hubbard U obtained from the SC 2-atom cell, 
once inserted in the 64-atom supercell, captures most of 
the effective interaction; second nearest neighbors shell 
brings some significant corrections to the final extrap- 
olated result, while third nearest neighbor shell has a 
smaller effect. We believe that contributions from fur- 
ther neighbor rapidly vanish and that an accurate value 
of U can be extracted from the SC supercell containing 
16 atoms. The extrapolation from this cell to larger cells 
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FIG. 3: Lattice spacing dependence of the calculated Hubbard 
U parameter for Iron. 

brings only minor variations which are within the finite 
numerical accuracy that we estimate within a fraction of 
an eV. From this analysis our estimate for the Hubbard 
U in elemental Iron at the experimental lattice parameter 
is therefore 2.2 ± 0.2 eV. 

This results is in very good agreement with the exper- 
imental estimates 0, ydl , but disagrees with Anisimov 
and Gunnarsson result [8j. We can only recall here that 
many technical details differ in the two approaches. In 
particular i) in the original approach the perturbed atom 
is disconnected from the rest of the crystal by removing 
all hopping terms, thus leaving a rather unphysical envi- 
ronment to perform the screening, while in our approach 
the actual system is allowed to screen the perturbation, 
ii) the Atomic Sphere Approximation (ASA) was em- 
ployed in the original LMTO calculation while no shape 
approximation is made in our case. 

In order to further test our approach on this element we 
investigate the dependence of the Hubbard parameter on 
crystal structure. The dependence of the calculated in- 
teraction parameter on the lattice spacing of the unit cell 
is shown in fig.|2lwhere a marked increase of the Hubbard 
U can be observed when the lattice parameter is squeezed 
below its experimental value. Despite this may appear 
counterintuitive, as correlation effects are expected to be- 
come less important when atoms gets closer, one should 
actually compare the increasing value of U with the much 
steeper increase of bandwidth when reducing the inter- 
atomic distance. Upon increase of the lattice parameter 
the Hubbard parameter should approach the atomic limit 
that can be estimated from all-electron atomic calcula- 
tions where the local neutrality of the metallic system is 
maintained: U = E{d 8 s°) + E{d e s 2 ) - 2 x E(d 7 s 1 ) = 2.1 
eV, in reasonable agreement with the results of fig. |21 
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TABLE I: Comparison between the calculated lattice con- 
stant (ao), bulk modulus (-Bo) and magnetic moment (/io) 
within several approximate DFT schemes and experimental 
results quoted from [28|. 



q (a.u.) Bp (Mbar) fxp (/xg) 
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1.68 
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2.33 


2.10 
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2.12 


2.60 


5.34 


1.53 


2.00 
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Using the calculated volume dependent Hubbard U pa- 
rameter we have studied the effect of the LDA+U ap- 
proximation on the structural properties of Iron. Results 
are reported in table [I] where they are compared with re- 
sults obtained within LSDA and ct-GGA(PBE) approxi- 
mation and with experimental data. From these data it 
appears that, although simple er-GGA(PBE) approxima- 
tion appears to be superior in this case, LDA+U provides 
a reasonable description of the data, of the same qual- 
ity as LSDA. In weakly correlated metals it has been 
suggested 29} that a formulation of LDA+U in terms 
of occupancy fluctuations around the uniform occupancy 
of the localized level could be more appropriate than the 
standard one. This "around mean field" (AMF) LDA+U 
approach has been revisited recently |3iJ,|3j| and an "op- 
timally mixed" scheme has also been proposed We 
don't want to enter in this discussion here, but we men- 
tion that by following the AMF recipe the description of 
structural and magnetic properties of metallic Iron im- 
proves as it is evident from tabled 

Using the calculated value of U we have obtained the 
electronic structure of Iron at the experimental lattice 
spacing. The theoretical band structure obtained using 
the AMF version of LDA+U is reported in fig. ^together 
with some experimental results |32j. The overall agree- 
ment is rather good for this scheme. However, when using 
the standard LDA+U scheme a somehow worse agree- 
ment with experimental data was obtained, mainly due 
to a rigid downward shift of the majority spin bands of 
about 1 eV. This is an indication that LDA+U approx- 
imation may still require some fine tuning in order to 
describe accurately both strongly and weakly correlated 
systems puj . 

Let us proceed to examine the Cerium case. Elemen- 
tal cerium presents a very interesting phase diagram with 
a peculiar isostructural a — 7 phase transition between 
a low volume (a) and a high volume (7) phase, both 
FCC. This phase transition has attracted much experi- 
mental and theoretical interest and in the last 20 years 
|33|. many interpretations have been put forward to ex- 
plain its occurrence. It is clear now that standard LDA 
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FIG. 4: Band structure of bulk iron obtained within the AMF 
LDA+U approach. Green lines are for minority spin states, 
black ones for majority spin levels. Photoemission results 
from |32| are also reported for comparison. 



or GGA approximations do not describe the transition 
and it appears that a treatment of the correlation at the 
DMFT level might be required j^H, however a full un- 
derstanding of the nature of the transition is still under 
debate |35|. Here, we do not want to address this deli- 
cate topics but we simply want to follow Anisimov and 
Gunnarsson 8] by computing the Hubbard U parameter 
for elemental cerium in the high volume 7 phase. 

The interaction of valence-electrons with Ce nuclei and 
its core electrons was described by a non-local ultrasoft 
pseudopotential 0] generated in the 5s 2 5p 6 5d 1 4f 1 elec- 
tronic configuration. Kinetic cutoffs of 30 Ry and 240 Ry 
were adopted for wavefunction and charge density Fourier 
expansion. The LSDA approximation was adopted for 
the exchange and correlation functional. Brihouin Zone 
integrations where performed using 8x8x8 Monkhorst 
and Pack special point grids |26l| using Methfessel and 
Paxton smearing technique j27j with a smearing width 
of 0.05 Ry. 

To obtain the response to an isolated perturbation we 
have perturbed a Cerium atom in three different cells: i) 
the fundamental face-centered cubic (FCC) cell contain- 
ing just one inequivalent atom, ii) a simple-cubic (SC) 
cell containing 4 atoms (giving access to the first nearest- 
neighbor response) and Hi) a 2x2x2 FCC cell (8 in- 
equivalent atoms) including also the response of second- 
nearest neighbor atoms. The result of these calculations 
and their extrapolation to very large SC cells is reported 
in Fig. where it can be seen that the converged value 
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FIG. 5: Calculated Hubbard U in metallic Cerium for 
different supercells. Lines connect results from the cell- 
extrapolation procedure and different symbols correspond to 
inclusion of screening contributions up to the indicated shell 
of neighbors of the perturbed atom. 




for U approaches 4.5 eV. 

The screening in metallic cerium is extremely localized, 
as can be seen from the fact that inclusion of the first- 
nearest neighbor response is all is needed to reach con- 
verged results. This is at variance with what we found in 
metallic Iron where third nearest-neighbor response was 
still significant (see Fig. (2J. The calculated value is not 
far from the value (5-7 eV) expected from empirical and 
experimental estimates |2fj, |23|, |24j , especially if we con- 
sider that the parameter U we compute plays the role 
of U — J in the simplified rotational invariant LDA+U 
scheme adopted [l4^ . 

As a check, we performed all-electron atomic calcula- 
tions for Ce + ions where localized 4/ electrons were pro- 
moted to more delocalized 6s or 5d states and obtained 
U = E(f 3 s°) + Eifs 2 ) - 2 x £(/V) = 4.4 eV, or 
U = E(ps°d 1 ) + Eifs 2 ^) - 2 x ^(/Vd 1 ) = 6.4 eV, 
depending on the selected atomic configurations. This 
confirms the correct order of magnitude of our calculated 
value in the metal. 

The present formulation is therefore able to provide 
reasonable values for the on-site Coulomb parameter 
both in Iron and Cerium, at variance with the original 
scheme of ref. where only the latter was satisfactorily 
described. We believe that a proper description of the 
interatomic screening, rather unphysical in the original 
scheme where atoms were artificially disconnected from 
the environment, is important to obtain a correct value 
for Hubbard U parameter, especially in Iron where this 
response is more long-ranged. 



FIG. 6: The unit cell of FeO: blue spheres represent Oxygen 
ions, red ones are Fe ions, with arrows showing the orienta- 
tion of their magnetic moments. Ferromagnetic (111) planes 
of iron ions alternate with opposite spins producing type II 
antiferromagnetic order and rhombohedral symmetry. 



Transition metal monoxides: FeO and NiO 

The use of the LDA+U method for studying FeO is 
mainly motivated by the attempt to reproduce the ob- 
served insulating behavior. In fact, as for other tran- 
sition metal oxides (TMO), standard DFT methods, as 
LDA or GGA, produce an unphysical metallic character 
due to the fact that crystal field and electronic struc- 
ture effects are not sufficient in this case to open a gap 
in the three-fold minority-spin ti g levels that host one 
electron per Fe 2+ atom. As already addressed in quite 
abundant literature on TMO (and FeO in particular), a 
better description of the electronic correlations is neces- 
sary to obtain the observed insulating behavior and the 
structural prop erties of this compound at low pressure 
HE l3Sll39j. The application of our approach to this 
material will thus allow us to check its validity by com- 
parison of our results with the ones from experiments and 
other theoretical works. 

The unit cell of this compound is of rock-salt type, 
with a rhombohedral symmetry introduced by a type II 
antiferromagnetic (AF) order (see fig. |f)J) which sets in 
along the [111] direction below a Neel temperature of 198 
K, at ambient pressure. 

The calculations on this materials were all performed 
in the antiferromagnetic phase starting from the cubic 
(undistorted) unit cell of fig. with the experimental 
lattice spacing. We used a 40 Ry energy cut off for the 



10 



5.0 



4.0 



3.0 



2.0 



e — oci 

□ □ C4 

OC16 



0.0 50.0 100.0 150.0 200.0 

Number of Fe ions per unit cell 



250.0 



> 

m 

>^ 
tot) 

<d 

a 
w 




FIG. 7: Convergence of Hubbard U parameter of FeO with 
the number of iron included in the supercell used in the ex- 
trapolation. Lines connect results including the screening con- 
tributions extracted from the indicated cell. 



electronic wavefunctions (400 Ry for the charge density 
due to the use of ultrasoft pseudopotentials [13( both for 
Fe and O) and a small smearing width of 0.005 Ry which 
required a 4x4x4 k-points mesh. 

To compute the Hubbard effective interactions, we per- 
formed GGA calculations with potential shifts on one 
Hubbard site in larger and larger unit cells, that we 
named CI, C4, and C16, containing 2, 8, and 32 iron 
ions respectively, and extrapolated their results up to 
a supercell containing 256 magnetic ions (called C128). 
The result for the undistorted cubic cell at the exper- 
imental lattice spacing is reported in fig. [7| We can 
observe that the effective interaction obtained from C4 
is already very well converged, when extrapolated to the 
largest cell, with respect to inclusion of screening from 
additional shells of neighborers, 

The final result for the Hubbard U is 4.3 eV which is 
smaller than most of the values obtained (or simply as- 
sumed) in other works [23, HE Ei3 • If we use this value 
in a LDA+U calculation we can obtain the observed in- 
sulating behavior as shown in the band structure plot of 
fig. |H1 where a comparison is made with GGA (metallic) 
results. 

A gap opens around the Fermi level whose minimal 
width is about 2 eV. The band gap is direct and located 
at the r point. The corresponding transition, of 3d(Fe)- 
2p(0)— >4s(Fe) character, should be quite weak due to 
the vanishing weight of Iron s states at the bottom of 
the valence band (fig. bottom picture). We can ex- 
pect that a stronger absorption line will appear instead 
around 2.6 eV due to the transition, of 3d(Fe)-2p(0) — » 
3d(Fe) character, among two pronounced peaks of the 
density of states around the Fermi level. This picture 
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FIG. 8: The band structure of FeO in the undistorted (cubic) 
AF configuration at the experimental lattice spacing obtained 
within GGA (top panel) and LDA+U using the computed 
Hubbard U of 4.3 eV (bottom panel). The zero of the energy 
is set at the top of the valence band. 



is in very good agreement with experiments (and other 
theoretical results |3tA |40j ) where a first weak absorp- 
tion is reported between 0.5 and 2 eV and a stronger line 
appears around 2.4 eV (il). The large mixing between 
majority-spin Iron 3d states and the Oxygen 2p manifold 
over a wide region of energy and the finite contribution 
of the Oxygen states at the top of the valence band — a 
feature not present within ct-GGA (see top panel in fig. 

— are also in good agreement with experiments, which 
indicate for FeO a moderate charge transfer character of 
the insulating state. 

Despite our U is smaller than the ones used in litera- 
ture, we find a good agreement of our results about the 
electronic structure of the system with experiments and 
other theoretical works. These findings confirm the va- 
lidity of our internally consistent method to compute U. 
We now want to extend its application to the study of 
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FIG. 9: Projected density of states of FeO in the undis- 
torted (cubic) AF configuration at the experimental lattice 
spacing obtained within GGA (top panel) and LDA+U using 
the computed Hubbard U of 4.3 eV (bottom panel). 



structural properties. This is indeed a very important 
test because a good ab-initio method should be able to 
describe the true ground state of a system and provide 
a complete description of both electronic and structural 
properties. Furthermore the plane-wave implementation 
we use allows a straightforward calculation of Hellmann- 
Feynman forces and stresses, thus giving easily access to 
equilibrium crystal structure. 

As observed in experiments [42(, the cubic rock salt 
structure of FeO shown in fig. becomes unstable un- 
der a pressure of 16 GPa (at room temperature) toward 
a rhombohedral distortion. In the distorted phase the 
unit cell is elongated along the [111] direction with a 
consequent shrinking of the interionic distances on the 
(111) planes. This transition is driven by the onset of the 
AFII magnetic order (the Neel temperature reaches 
room value at about 16 GPa) which imposes a rhombohe- 
dral symmetry even in the cubic phase. Upon increasing 
pressure above the threshold value the distortion of the 
unit cell is observed to increase producing more elongated 
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FIG. 10: The pressure dependence of the rhombohedral angle 
in FeO for the various approximations described in the text is 
compared with experimental results. These latter results were 
extracted extrapolating the data for the non stoichiometric 
compound Fei-xO up to the stoichiometric composition [I^, 
0. 



structures |42(. 

We have computed the Hubbard U on a grid of possible 
values for the rhombohedral distortion and cell parame- 
ter and then from the corresponding total energy calcu- 
lations we determined the rhombohedral distortion and 
the enthalpy of the system as a function of the pressure 
up to 250 Kbar. 

As evident from fig. ^3 while GGA overestimates the 
rhombohedral distortion and his pressure dependence, 
LDA+U method — in the standard electronic configura- 
tion examined so far-overcorrects the GGA results and 
introduces even larger errors with respect to experimen- 
tal results. In fact not only we obtain a distortion with 
the wrong sign (of compressive character along the [111] 
direction), but also the wrong pressure dependence. The 
reason for this failure can be traced back to the different 
occupation of the orbitals around the gap/Fermi level in 
the two cases. Even in the undistorted cell, the rhombo- 
hedral symmetry, induced by the antiferromagnetic or- 
der, lifts the degeneracy of the minority spin t2 g states 
of iron and split them in one state of A\ g character — 
which is essentially the m=0 (z 2 ) state along the [111] 
quantization axis — and two states of e g symmetry local- 
ized on the iron (111) planes. Within GGA, the Iron 
minority-spin 3d electrons partially occupy the two equiv- 
alent e g orbitals giving rise to two half filled bands and a 
(wrong) metallic state which is delocalized on the (111) 
plane. The system gains energy by filling the lowest half 
of the e g states and tends to elongate in the [111] direc- 
tion, shrinking in the plane, because this increases the 
overlap of the e g states and their bandwidth. Within 
LDA+U, fractional occupation of orbitals is energetically 
disfavored and the system would like to have completely 
filled or empty 3d states. In the standard unit-cell con- 
sidered so far in the literature — and used by us in the 
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FIG. 11: Lattice distortion in the (111) iron planes used to 
induce symmetry breaking in the electronic configuration of 
FeO. 



calculation above — this can be accomplished only by fill- 
ing the non-degenerate A\ g level, corresponding to wave- 
functions elongated along [111], and pushing upward in 
energy the in-plane e g states, leaving them empty. As 
a consequence, the system tends to pull apart the ions 
on the same (111) plane, so that the bandwidth of the 
state in the plane is reduced, and increases instead the 
inter-plane overlap of the A\ g states. This simple picture 
gives an explanation of the fact that GGA overestimates 
the elongation of the unit cell in the [111] direction, as 
well as the (wrong) compressive behavior of the standard 
LDA+U solution. We are thus left with the paradoxical 
situation that a correct pressure dependence of the struc- 
tural properties can be obtained from the wrong band 
structure and viceversa. 

We have found that it is possible to solve this para- 
dox by allowing the possibility that the system partially 
occupies, as within GGA, the e g levels, thus maintaining 
the driving force for the right rhombohedral deformation, 
and still opens a gap, as in standard LDA+U, by some 
orbital ordering that breaks the equivalence of the iron 
ions in the (111) plane. This possibility has been some- 
times proposed in literature |39l but has never been 
clearly addressed. 

From a simple tight-binding picture one finds that the 
optimal broken symmetry phase would be the one where 
occupied e g orbitals have the highest possible hopping 
term with unoccupied e g orbitals in nearest-neighbor 
atoms in the plane, in order to maximize the kinetic 
energy gain coming from derealization, and the low- 
est possible hopping term with neighboring occupied e g 
orbitals, in order to minimize bandwidth that tends to 
destroy the insulating state. In bipartite lattice this is 
simply achieved by making occupied orbitals in nearest- 
neighbor sites orthogonal but, in the triangular lattice, 
formed by iron atoms in (111) planes, this is not exactly 
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FIG. 12: The projected density of states of FeO as obtained 
in the "standard" LDA+U ground state (top panel) and in 
the proposed broken symmetry phase (bottom panel) . On the 
right of each DOS is a picture of the corresponding occupied 
Fe-3d minority states. 



possible, the system is topologically frustrated and some 
compromise is necessary. 

It is generally believed ji^ that Heisenberg model in 
the triangular lattice, to which our system resemble in 
some sense, displays a three-sublattice 120° Neel long- 
range order. We thus imposed a symmetry breaking to 
the system where three nearest-neighbor atoms in the 
(111) plane were made inequivalent by slightly displac- 
ing them from the ideal positions in the way shown in 
fig. 1111 This induced the desired symmetry breaking of 
the electronic structure and opened a gap that was ro- 
bust and persisted when the atoms were brought back 
into the ideal positions. We found, quite satisfactorily, 
that the new broken symmetry phase (BSP) corresponds 
to a lower energy minimum than the " standard" LDA+U 
solution and that therefore it is, to say the least, a more 
consistent description of the ground state of FeO. The one 
depicted in fig. ^]is, of course, only one of three equiva- 
lent distortions we could have imposed to the electronic 
structure of the system and three symmetry related BSPs 
could be defined. In the actual system an effective equiv- 
alence of the ions in the (111) planes is probably restored 
by a (dynamical) switching among equivalent states but 
considering the atoms as strictly equivalent, as in the 
standard solution, leads to incorrect results. 

The comparison of the projected density of state in the 
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" standard" LDA+U solution and in the novel BSP phase 
is shown in fig. 1121 where also a pictorial representation 
of the occupied minority-spin orbitals in the two cases 
is shown. As we can observed, no remarkable qualita- 
tive difference in the DOS appears apart from the differ- 
ent ordering of the d states around the gap. In fact the 
minority-spin d electron is now accommodated on a state 
lying on the (111) plane (shown on the right panel) while 
the one with A\ g (z 2 ) character has been pushed above 
the energy gap. The gap width and the charge transfer 
character of the system do not change significantly and 
are still in very good agreement with the experiments. 

We repeated the structural calculations (according to 
the same procedure described above) in the BSP, and 
obtained the LDA+U (BSP) curve reported in fig. HOI 
The agreement with experiments is much improved with 
respect to both GGA and LDA+U "standard" ground 
states. The mechanism leading to the pressure behavior 
in BSP case is basically the same already producing the 
correct evolution of distortion in the GGA calculations. 
When the unit cell elongates along the cubic diagonal 
the iron ions in the (111) plane get closer and the hop- 
ping between nearest-neighbor orbitals increased with a 
consequent lowering of the electronic kinetic energy. 

We therefore conclude that LDA+U, not only improves 
the description of the structural and electronic properties 
with respect to GGA, but that a close examination of 
both electronic and structural properties is in this case 
necessary in order to describe the correct ground state of 
the system. 

Another classical example of TMO we want to study 
in order to test the present implementation of LDA+U 
is Nickel Oxide. It is a very well studied material and 
there is a good number of theoretical 01 an d experi- 
mental works, including some photoemission experiments 
H^H^, our results can be compared with. At variance 
with FeO, no compositional instability is observed for 
NiO so that the stoichiometric compound is easy to study 
and is much better characterized than iron oxide. It has 
cubic structure with the same AF spin arrangements of 
rhombohedral symmetry as FeO, but does not show ten- 
dencies toward geometrical distortions of any kind and is 
therefore easier to study. 

In this case we did not perform any structural relax- 
ation and calculated the value of U at the experimental 
lattice spacing for the cubic unit cell imposing the rhom- 
bohedral AF magnetic order which is the ground state 
spin arrangement for this compound. The GGA approx- 
imation (in the PBE prescription) was used in the calcu- 
lation. US pseudopotentials for Nickel and Oxygen (the 
same as in FeO) were used with the same energy cutoffs 
(of 40 and 400 Ry respectively) for both the electronic 
wavefunctions and the charge density as for FeO and also 
the same 4x4x4 k-point grid for reciprocal space integra- 
tions. 

In the calculation of the Hubbard U of NiO we did not 



studied the convergence properties of U with system size 
as we did in FeO but, assuming a similar convergence 
also in this case, we performed a constrained calculation 
only in the C4 cell and then extrapolated the obtained 
result to the C128 supercell. The calculated value of the 
U parameter is 4.6 eV. This value is smaller than litera- 
ture values for the same parameter that are rather in the 
ran ge o f 7-8 eV , however it has been recently pointed 
out [13, LLJ] that in obtaining these values self-screening 
of d electrons is neglected and that better agreement with 
experimental results is obtained using an effective Hub- 
bard U of the order of 5-6 eV. 

The magnetic moment of the Ni ions is correctly de- 
scribed within the present GGA+U approach which gives 
a value of 1.7 fis well within the experimental range of 
values ranging from 1.64 and 1.9 \ib |48l l49l| . better than 
the value of 1.55 obtained within GGA. 

In fig. El and fig. E] the band structure and atomic- 
state projected density of states of NiO obtained with this 
value of U is shown, along with the results of standard 
GGA, and compared with the photoemission data in the 
rX direction extracted from ref . [46l l47j | . 

Despite the agreement with the experimental band- 
dispersion is not excellent — the valence band width is 
somehow overestimated by both GGA and GGA+U 
calculations — , GGA+U band structure reproduces well 
some features of the photoemmission spectrum for this 
compound and gives a much larger band gap than the 
one obtained within GGA approximation. A very im- 
portant feature to be noticed in the density of states re- 
ported in fig. El is the fact that GGA+U modifies quali- 
tatively the nature of the states at the top of the valence 
band, and hence the nature of the band gap: in GGA 
approximation the top of valence band is dominated by 
Nickel cZ-states while in the GGA+U calculation the Oxy- 
gen p-states give the most important contribution. In 
both approaches the bottom of the conduction band is 
mainly Nickel d-like and therefore the predicted band 
gap is primarily of charge-transfer type within GGA+U, 
in a gree ment with experimental and theoretical evidence 
|ifl l5fi Iftlj . while it is wrongly described as of Mott- 
Hubbard type according GGA approximation. 

Our GGA+U value for the optical gap is « 2.7 eV 
around the T point, smaller that commonly accepted 
experimental values that range from 3.7 to 4.3 eV 
5 2.l53ll54ll55| . More recently however, a re-examination 
56| of the best available optical absorption data 
pointed out that optical absorption in NiO starts at pho- 
ton energy as low as 3.1 eV, not far from our theoretical 
result. Indeed, Bengone and coworkers reported re- 
cently an LDA+U calculation in NiO where different em- 
pirical values of U were employed. When U — 5 eV was 
used — a value close to our present first-principles result — 
, they obtained an optical gap of 2.8 eV, very close to 
our results, and an excellent agreement between the cal- 
culated and experimental [52^ optical absorption spectra. 
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FIG. 13: The band structure of NiO in the undistorted (cu- 
bic) AF configuration at the experimental lattice spacing ob- 
tained within GGA (top panel) and with the computed Hub- 
bard U of 4.6 eV (bottom panel). The zero of the energy is 
set at the top of the valence band. Experimental data from 
ref. |4rt (empty symbols) and [i^ (solid symbols) are also 
reported. 



The same calculation with the literature value of U = 8 
eV, gave a larger value for the optical gap but a very poor 
agreement with the experimental absorption spectrum. 



Minerals: Fayalite 

As a final example we want to apply the present 
methodology to Fayalite, the iron-rich end member of 
(Mg,Fe) 2 Si04 olivine (orthorhombic structure), one of 
the most abundant minerals in Earth's upper mantle. 
Recently [5^ we showed that, although good structural 
and magnetic properties could be obtained for this min- 
eral within LDA or GGA, its electronic properties were 
incorrectly described as metallic, confirming the corre- 
lated origin of the observed insulating behavior. 




LDA+U 
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FIG. 14: Projected density of states of NiO in the undistorted 
AF configuration at the experimental lattice spacing obtained 
with U = 4.6 eV. 



From x-rays diffraction studies it is known that Fay- 
alite has an orthorhombic cell, whose experimental lattice 
parameters are (in atomic units) a = 19.79, b = 11.50, 
c = 9.11. The unit cell (depicted in fig. 1151) contains 
four formula units, 28 atoms: 8 iron, 4 silicon, and 16 
oxygen atoms. Silicon ions are tetrahedrally coordinated 
to oxygens, whereas iron ions occupy the centers of dis- 
torted oxygen octahedra. The point group symmetry of 
the non magnetic crystal is mmm (D211 in the Schoenflics 
notation) and the space group is Pnma. The magnetiza- 
tion of iron reduces the original symmetry and only half 
of the symmetry operations survive. The general expres- 
sion for the internal structural degrees of freedom is given 
in tabic [D] in the Wyckoff notation [Hi} . 

Iron sites can be divided into two classes (see fig. IT51 
and tab.nj: Fel centers which are structured in chains 
running parallel to the b, [010], side of the orthorhom- 
bic cell, and Fe2 sites which belong to mirror planes for 
the non magnetic crystal structure perpendicular to the 
b side and cutting it at 1/4 and 3/4 of its length. The 
main structural units are the iron centered oxygen octa- 
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FIG. 15: The unit cell of Fayalite. Large dark ions are Fe, 
small dark ions are O, light ions are Si. 

TABLE II: Definition of the Wyckoff structural parameters 
appropriate for Fayalite structure 

Ion Class Coordinates 

Fel 4a (0,0,0), (1/2,0,1/2) 

(0,1/2,0), (1/2,1/2,1/2) 
Fe2, Si, 01, 02 4c ±(u,l/4,v), 

±(u+l/2,l/4,l/2-v) 
03 8d ±(x,y,z), ±(x,l/2-y,z), 

±(x+l/2,l/2-y,l/2-z), 
±(x+l/2,y,l/2-z) 



hedra which are distorted from the cubic symmetry and 
tilted with respect to each other both along the chains 
and on nearest Fe2 sites. Fayalite is known to be an anti- 
ferromagnetic (AF) compound with slightly non collincar 
arrangement of spin on Fel iron site (this non collincar- 
ity will not be addressed here) . Magnetic moments along 
the central and the edge Fel chains are antiferromagnet- 
ically oriented and from our previous work |57| the most 
stable spin configuration is the one in which the magne- 
tization of Fe2 ion is parallel to the one of the closest 
Fel iron. This magnetic structure is consistent with an 
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FIG. 16: Convergence of Hubbard parameters of Fayalite with 
the number of iron included in the supercell used in the ex- 
trapolation. Ul is the value obtained for Fel ions, U2 the one 
for Fe2. 



iron-iron magnetic interaction via a superexchange mech- 
anism through oxygen p orbitals. 

The calculation of U was performed for the experimen- 
tal geometry, in the above mentioned spin configuration. 
As the primitive unit cell of fayalite is already quite large, 
we performed the constrained calculation only in this cell 
and used larger supercells only to extrapolate the results. 
We considered three supercells in addition to the prim- 
itive one: i) a cell duplicated in the [0, 1, 0] chain direc- 
tion (a 1x2x1 supercell), containing 16 iron atoms; ii) 
a cell, containing 64 iron ions, obtained by duplicating 
the primitive structure in all directions (a 2x2x2 super- 
cell) and in) a 4x4x2 supercell (256 iron ions). Other 
computational details where similar to those used in our 
previous work |57| . As GGA approximation provided a 
slightly better description of the system than LDA, we 
assumed this functional as the starting point to be im- 
proved; the same pseudopotentials used in ref. jH^I for 
Fe, O and Si were adopted here; somehow larger energy 
cutoff for the electronic wave functions and charge den- 
sity (36 and 288 Ry respectively) and a small smearing 
width of 0.005 Ry were used. A 2x4x4 Monkhorst-Pack 
grid of k-points in the primitive cell was found sufficient 
for the BZ integration. 

The results of the U calculation for the two different 
families of iron sites (Fel and Fe2) are reported in fig. 
1161 where the rapid convergence with respect supercell 
dimension can be seen. The final results for the on-site 
Coulomb parameters are U\ = 4.9 eV for Fel ions and 
U2 = 4.6 eV for Fe2, which are in fairly good agreement 
with the approximate (average) value of 4.5 eV obtained 
in ref. |57| from a rather crude estimate. 

The GGA+U band structure of Fayalite is shown in fig. 
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FIG. 18: Some atomic-projected density of states of Fayalite 
obtained within the present LDA+U approach. Contributions 
from majority- and minority-spin 3d states of one of the Fel 
iron ions and from the total 2p manifold of one oxygen ion 
are shown. 



FIG. 17: The band structure of Fayalite obtained within the 
present LDA+U approach. The zero of the energy is set to 
the top of the valence band. Complete degeneracy among 
spin up and spin down states is present. 



1 171 while in fig.^3]some atomic-projected density of states 
are reported. At variance with the GGA results reported 
in ref. [H3 a band gap of about 3 eV now separates the 
valence manifold from the conduction one, in reasonable 
agreement with the experimental result of about 2 eV 
|59j at zero pressure. 

The minority spin t2 g manifold of iron ions, that within 
GGA cross the Fermi energy, is split into two subgroups 
by the gap opening. The conduction-band states are 
shrunk to a narrow energy range and moved above the 
bottom of the iron s-states band which remains almost 
unaffected; the lower-energy minority-spin ri-states, in- 
stead, merge in the group of states below the Fermi level 
where they mix strongly with states originating from 
Oxygen p orbitals: the two sets of states, well separated 
in the GGA results, collapse into a unique block. The 
most evident consequence of the gap opening consists in 
a pronounced shrinking of the d states of iron which be- 
come flatter than in the GGA case. This is evident on 
the top of the valence band, but also for states well below 
this energy level, which thus reveal a more pronounced 
atomic-like behavior. Beside the gap opening between 
the two groups of the minority-spin states, a strong mix- 
ing occurs among the oxygen p states and the iron d 
levels over a rather large region extending down to 8 eV 



TABLE III: Comparison of the experimental and LDA+U 
calculated values for the Wyckoff structural parameters of 
Fayalite as defined in table ITT1 



Ion 


u 


V 


X 


y 


z 


Exp. 
















Fe2 


0.780 


0.515 










Si 


0.598 


0.071 










01 


0.593 


0.731 










02 


0.953 


0.292 










03 






0.164 


0.038 


0.289 


GGA + U 
















Fe2 


0.779 


0.515 










Si 


0.597 


0.072 










01 


0.593 


0.735 










02 


0.951 


0.289 










03 






0.165 


0.036 


0.286 



below the top of the valence band. A finite contribu- 
tion of the oxygen states is present close to the top of 
the valence manifold showing that the gap is mainly of 
Mott-Hubbard type with a partial charge-transfer char- 
acter. 

We have then relaxed the geometric structure of the 
system (both internal and cell degrees of freedom) assum- 
ing no dependence of 17% and P2 on the atomic configu- 
ration. The resulting structural parameters (a — 20.18, 
b = 11. 75, c — 9.29 atomic unit) as well as the internal co- 
ordinates reported in table llll| are in very good agreement 
with the experimental results, even better than the al- 
ready satisfactory agreement obtained in ref. JH^I within 
GGA. 
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Although wc did not studied other spin-configurations, 
magnetic properties seam to improve slightly in the 
GGA+U approximation. The magnetic moment on each 
iron (both Fel and Fe2) was found to be 3.9 in closer 
agreement with the spin-only value (4 /xs) of the exper- 
imental result (4.4 tis) than the one obtained by GGA 
only (3.8 /is). This improvement is probably due to the 
enhanced atomic-like character of iron d states, which is 
consequence of the gap opening. 

In conclusion, the GGA+U provides a quite good de- 
scription of structural, magnetic and electronic properties 
of fayalite, reproducing the observed insulating behavior 
with a reasonable value for its fundamental band gap. 

SUMMARY 

In this work we have reexamined the LDA+U approxi- 
mation to DFT and a simplified rotational-invariant form 
of the functional was adopted. We then developed a 
method, based on a linear response approach, to calculate 
in an internally consistent way the interaction param- 
eters entering the LDA+U functional, without making 
aprioristic assumption about screening and/or basis set 
employed in the calculation. Our methodology was then 
successfully tested on a few systems representative of nor- 
mal and correlated metals, simple transition metal oxides 
and iron silicates. In all cases we obtained rather accu- 
rate results indicating that our scheme allows us to study 
both electronic and structural properties of strongly cor- 
related material on equal footing, without resorting to 
any empirical parameter adjustment. 

This work has been supported by the MIUR under the 
PRIN program and by the INFM in the framework of the 
Iniziativa Trasversale Calcolo Parallelo. 
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